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We derive an alternative representation for the /-electron spectral function of the Falicov-Kimball model from 
the original one proposed by Brandt and Urbanek. In the new representation, all calculations are restricted 
to the real time axis, allowing us to go to arbitrarily low temperatures. The general formula for the retarded 
Green's function involves two determinants of continuous matrix operators that have the Toeplitz form. 
By employing the Wiener-Hopf sum equation approach and Szego's theorem, we can derive exact analytic 
formulas for the large-time limits of the Green's function; we illustrate this for cases when the logarithm of 
characteristic function (which defines the continuous Toeplitz matrix) does and does not wind around the 
origin. We show how accurate these asymptotic formulas are to the exact solutions found from extrapolating 
matrix calculations to the zero discretization size limit. 
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1. Introduction 

The Falicov-Kimball model [ [T] has often been studied as one of the simplest models of strong 
electron correlations. Indeed, the exact solution in the limit of large dimensions displays charge 
density wave order, the Mott metal-insulator transition, and phase separation [ 2J. It is particularly 
relevant to present a paper on this model in an issue dedicated to Prof. Stasyuk, as he has worked 
on this model and closely related ones throughout his career; this began with his first article on the 
Falicov-Kimball model [ 3 while his most relevant work to this contribution involves approximate 
solutions for the /-electron spectral functions [ [H [51 [5] . 

Our focus in this work is on the spectral function of the /-electrons in the exact solution of 
Falicov-Kimball model with dynamical mean-field theory This problem was first investigated by 
Brandt and Urbanek [ [7] and by Janis [ [8] . The numerical approach of Brandt and Urbanek was 
extended by Zlatic et al. [ f5j and Freericks, Turkowski and Zlatic [ [TOl HTJ |T^ within the origi- 
nal Brandt-Urbanek framework. An alternative approach was taken by employing the numerical 
renormalization group by Anders and CzychoU [ [T3| . Some analytic work was completed by Liu [ 
114) which was related to early work on this problem [ I15j and proceeded from the perspective of 
the X-ray edge problem [ I16|, I17j . Here, we develop a new approach which is numerically more 
tractable than other approaches and allows us to develop asymptotically exact formulas for the 
Green's function in the time representation. 

The single-site Hamiltonian of the Falicov-Kimball model [ 1 involves two types of electrons — 
conduction electrons, denoted by the letter d and localized electrons, denoted by the letter / 
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Hioc = EfUf + UudUf - fj,{nd + Uf) = Ho{l - Uf) -\- HiUf 
Ho = -y.nd, 

Hi = Ef - /i + ([/ - ^i)nd, 



(1.1) 
(1.2) 
(1.3) 
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Figure 1. The Kadanoff-Baym-Keldysh contour. The contour starts at time 0, moves along 
the real time axis to time t, moves back along the real axis to time 0, then extends down the 
imaginary axis to time — i/3. 

where = d)d and Uf = f^f are the number operators for the d- and / electrons, respectively, U 
is the on-site Coulomb interaction between the d and / electrons, and Ef is the site-energy of the 
/ state. The full Hamiltonian on the lattice includes a repeat of this local Hamiltonian for each 
lattice site and a hopping of the conduction electrons between nearest-neighbor sites. The density 
matrix of the single-impurity problem is equal to 

p = e-f^"'°-%ei^^i^~i jdt' jdt''d\t')X^{t',t'')d{t'')^ 1, (1.4) 

where the time-ordering and integration are performed over the Kadanoff-Baym-Keldysh contour [ 
[TSl [T^ and (i — \ /T is the inverse temperature; the Kadanoff-Baym-Keldysh contour is shown 
in Fig. [T] — it starts at t = 0, runs out along the real axis to a maximal time tmax, then returns 
along the real axis back to i = 0, and finally runs along the negative imaginary axis down to —ip. 
Here, the dynamical mean- field Xc{t' , t") and chemical potential fj, are taken from the equilibrium 
solution of the conduction electron problem via dynamical mean-field theory [[2]. In other words, 
the dynamical mean field satisfies 

oo 

Ac(t,i') = -^ / dujlmX{u;)e-'^^'-''^[f{Lo)^e,it,t% (1.5) 

— oo 

where X{uj) is the ordinary dynamical mean field on the real axis, extracted from the dynamical 
mean-field theory solution of the model, f{uj) — 1/[1 -I- exp{(3uj)] is the Fermi-Dirac distribution 
function, and 8c(t, i') is the Heaviside function on the contour, equal to 1 when t is ahead of t' 
on the contour, equal to when t is behind t' on the contour and equal to 1/2 when t = t' on the 
contour (note that this definition of the Ac field is i times the definition in Ref. [H]). The partition 
function for the single-site problem contains two terms: 



Z = Tr e-'^J^iooTc exp \^~^ J^dt' dt" d\t')Xc{t' , t")d{t") 
which correspond to the different occupations of the /-particles, as follows: 



Zo + Zi, (1.6) 



Zo^[l-.e^iri(l--^)^ n,^0; (1.7) 

Z.^eP^^-^'fAl + e^^^-'^AWil ^ V nf = l. (1.8) 

Here, iujm — i'!TT(2m + 1) is the fermionic Matsubara frequency and 

P 

Xrn = /dre'"""A(r) (1.9) 
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is the dynamical mean field evaluated at the mth Matsubara frequency [the Ac field is a function 
only of the difference of the two time arguments when both times lie on the imaginary part of the 
Kadanoff-Baym-Keldysh contour and we use the notation A(t) = —iXc{—iT, 0), which corresponds 
to the conventional definition of the dynamical mean field for imaginary time] . 



2. Real-time Green's functions 

We consider the contour ordered Green's function for the /-electrons (with times t and t' both 
lying on the contour) 

G%t,t') = -i{TJ{t)p{t')), (2.1) 

where the time ordering is taken along the contour, the angle brackets denote taking a trace 
(weighted by the density matrix for the equilibrium system), and the time dependence of the 
fermionic operators is with respect to the local Hamiltonian -ffioc- Depending on the location of 
the time arguments on the contour, we obtain different real-time Green's functions: 
the greater Green's function 



G>{t - 1') = -i (fit) Pit')) = -4 E <p \f\ p) 



pq 



the lesser Green's function 

G<(t-i') = i(/nO/(i))=4E' 



{q\P\p){p\f\q) e'^'--'" 



i{ep—eg){t—t') . 



(2.2) 



(2.3) 



the time-ordered Green's function 

G}it - t') = -i {%f{t)f\t')) = e{t - t')G>{t - t') + e{t' - t)G<{t - t'); (2.4) 
and the antitime-ordered Green's function 

G}{t - t') = -i {Td{t)f\t')) = e{t' - t)G>{t - t') + Q{t - t')G<{t - t'). (2.5) 

Here \p) and \q) denote many-body eigenstates of the local Hamiltonian with eigenvalues Sp and 
£g, respectively, and Qit) is the ordinary Heaviside function. We are also interested in the retarded 
and advanced Green's functions, 



G}{t - 1') = -te{t - t'){[f{t), fHt')]+) = e(t - 1') [G>{t - 1') -G<{t- 1') 

G}{t - t') = iQ{t' - t) ( [fit), Pit')] \ = Bit' - t) \Gfit - t') -G>it- t') 



(2.6) 
(2.7) 

which are easily constructed from taking combinations of the lesser and greater Green's functions 
(note the symbol [...]-|- denotes the anticommutator of the two operators in the brackets). In 
equilibrium, the Green's functions depend only on the time difference of the real-time variables 
and the greater and lesser Green's functions also satisfy 



G>it~t') =-G>fit'-t), Gfit-t') =-Gfit'-t) 



(2.8) 



It is often useful to represent quantities in terms of frequencies, instead of timc\ The Fourier 
transforms of the greater and lesser Green's functions (to frequency space) are equal to: 



G^(w) = -27rz[l - /(w)]^/(w), GfiLo) = 27ri/(w)^/(w), 



where 



-fie, 



"jiplPq) {q\P\p)Siio~ 



(2.9) 
(2.10) 
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is the /-electron spectral density. The corresponding Fourier transforms of the retarded and time- 
ordered Green's functions are the following: 



uj — oj' + id ■' 



(2.11) 



and 



G'jiu;) ^ReG}{uj)+i Im G} (cj) tanh ^ . 



(2.12) 



We use the Kadanoff-Baym-Keldysh approach to calculate the real-time Green's function be- 
cause there is no simple way to perform the analytical continuation of the Matsubara frequency 
Green's function to the real axis [ [7] . We are interested in the retarded Green's function, so we 
consider the greater Green's function for positive time t > 



G>{t) 



~i^Tr<j e^'^^'^-^Tcexp 



i dt' dt"d\t')Xc{t',t")d{t") 



Pit)f{o)>; 



(2.13) 



(2.14) 



we use the fact that /^(0)/(t) can be replaced by /^(i)/(0) in the trace due to the time-translation 
invariance that we have in equilibrium. The greater and lesser Green's functions can be found for 
other time values by using the relations in equation (|2.8p . 

The first step in solving for the Green's function is to write equations of motion for the /- 
operators: 



and the lesser Green's function for negative time t < (or t = ~t > 0) 
Gf{t) ^ t^Ti- S^e-'^"'°-'%exp -i j dt' j dt" d\t')\^{t' ,t")d{t"] 



dm 

dt 
dfHt) 
dt 



^~i[Und{t)+Ef- ^i]f{t), 
^i[Undit) + Ef-^j.]fiF), 



(2.15) 
(2.16) 



and substitute their solutions into equations (|2.13p and (|2.14p . which yields the following expres- 
sions for the greater Green's function {t > 0) 

G>{t) ^-i^Trje-^^'-Tcexp ~i Jjt' J^dt"d\t')X^{t' ,t")d{t") 

/(0)/t(0) 



- i J dt'Ucit, t')nd{t') - i{Ef - fi)t 
and for the lesser Green's function {t — ~t > 0) 



G<{t) =j^Tr|e-'^^'-Teexp -i j dt' j dt" d\t')\{t' ,t")d{t") 



+ i / dt'Ucit, t')nd{t') + i{Ef - fi)t 



/t(0)/(0) 



where Uc(t,t') is a new time-dependent field, which is nonzero only on the upper branch of the 
Kadanoff-Baym-Keldysh contour and is equal to U for t' £ [0, t] and is zero otherwise. It is because 
this Uc field is not time-translation invariant that we need to use the Kadanoff-Baym-Keldysh 
formalism for the analytic continuation. Now we take into account the fact that the /-particle 
occupation number Uf = f^f is a conserved quantity of the Hamiltonian; for the greater Green's 
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function we have a projection onto states without /-particles (n/ =0) while for the lesser Green's 
function we have a projection onto states with one /-particle [nf = 1), which yields: 



G>{t) = -i_e-*(^/-'^)* 



X Tr <^ e 



T^exp 



i dt' dt"d^t')Xcit',t")d{t")-i / dt'Uc{t,t')ndit') 



and 



(2.17) 



(2.18) 



X Tr <^ e^(A'-t')nrf^ 



i dt' dt"d''{t')Xc{t',t")d{t"),+i / dt'Uc{t,t')nd{t') 



exp 



where the trace is now over the d-electrons only. 

We next expand the contour-ordered-exponent into a series that we can ultimately sum exactly. 
To begin, we need to recall Wick's theorem, where the expectation value of any number of pairs 
of fermionic operators can be expanded in terms of products of the two-operator contractions 



d{ti)d\t2) = iga{ti,t2), 
ga{ti,t2) = -t{%d{h)dHt2))^ = te-''-^''-'''> [/(e„)-ec(ti,t2)] 



(2.19) 



when the Hamiltonian is quadratic in the fermionic operators. Here we have eq = — /i for sites with 
Uf — {a = 0) and ei = U — for sites with Uf = 1 (a = 1). The traces in equations (j2.17p and 
(j2.18p can then be written in the following diagrammatic expression 



1 + e- 



exp ■ 



1 

+ 2 



(2.20) 



where the arrows denote the zero-order Green functions gait' , t") in equation (|2.19p and the wavy 
lines denote a generalized time-dependent hopping (A-ficld) 



\t{t' ,t") = \{t' ,t") + u,{tX)5c{t' ,t") 

for the greater Green's function and we have to replace 

UAt,t') ^ -C/e(t,i') 



(2.21) 



(2.22) 



for the lesser Green's function. 

Motivated by approaches that involve the integration over a coupling constant, we modify our 
time-dependent field by introducing a dependence on some new parameter that we denote by x 



UrXt,t')~.Ucit,t'\x), 

with the following limiting behavior 

C/,(i,t'|0) = 0, Ucit,t'\l)^Ucit,t'). 



(2.23) 



(2.24) 



Next, we take the derivative of the diagrammatic series in equation (|2.20p with respect to x and 
find: 



dxi'"} J 



= - I dt'^^^^^^^^G^';\t',t'\x), 
dx 



(2.25) 
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where we introduced a parameter-dependent Green's function G^^\t' ^t"\x), which is the solution 
of the foUowing Dyson equation: 

G'^;^\t' ,t"\x) =g^{t\t")+ f dh I dt2ga{t',hyXt{h,t2\x)Gl;;\t2,t"\x) (2.26) 



=9c.{t\t") + J dti J dhgo^it' M)K{h,t2)G\^\t2,t"\x) 

+ jdtigJt'M)U,{tM\x)du\ti^t"\^)- 

This expression holds for the greater Green's function {a — Q). For the lesser Green's function, one 
has to use the substitution in equation (j2.22p and set a = \. Finally, we find explicit expressions 
for the greater 

G> (t) = -iwo exp I -i{Ef -^i)t- J dxj^ dh ^E^Mll^c^^) ^ i (2.27) 



and lesser 



G<{t)^iwiexp^i{Ef~fi)t + dx j dh'^^^^^^^^^G^p {ti,h\x) ) , (2.28) 

Green functions, where 

^i'o = (1 - n/) = „ ^° , = (uf) = (2.29) 

are the average densities of sites without (wq) and with (wi) /-electrons. 
We can rewrite equation (|2.26p in the following two forms: 

G^^\t',t"\x)=g^mt',t"\x) + fdh f dt2g^u'':{t',h\x)X,{h,t2)G\;^\t2,t"\x), (2.30) 



G\^\t',t"\x) = G^{t',t") + J dhG^{t',h)Ucit,h\x)G\;^\h,t"\x). (2.31) 

The first form emphasizes the lambda field as an effective self-energy, while the second emphasizes 
the U field as the effective self-energy. In the top equation, we introduced the auxiliary Green's 
function 

gtjl^t',t"\x)^g^{t',t") + [ dhg^{t',h)U,{t,h\x)g^mh,t"\x) (2.32) 



as first defined by Brandt and Urbanek [Jj, while we introduced the bare time-ordered Green's 
function for the effective medium 

G„(i',t") =5a(i',i")+ / dh I dt2gc.{t'M)K{tut2)G^{t2^"), (2.33) 



in the bottom equation. The bare time-ordered Green's function for the effective medium can be 
determined directly on the real axis via 

Ga{t'X) j doj Im^„(a;)e-*'-(*'-*") [J{iv) - Q^t' ,t")] , (2.34) 

— OO 

with 

= — ^ w (2-35) 
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being the retarded effective medium on the real frequency axis. 

The first representation (|2.30p was originally used by Brandt and Urbanek [\7\ and improved by 
Zlatic et al. and Freericks, Turkowski and Zlatic [[TO] (see for review Freericks and Zlatic [[5]) 
and requires the calculation of continuous matrix operator determinants over the entire contour. 
In what follows, we shall instead use the second representation in equation ()2.3ip because in 
equations (I2.27P and (|2.31|) all times in Gij{ti,ti\x) are only on the real axis {0 < ti < t) [or 
(0 < ii < t)] and we can replace the contour integral in equation (|2.3ip by an integral over the 
upper real time branch of the Kadanoff-Baym-Keldysh contour where the time-dependent field 
Uc{t^ti\x) is nonzero 



-t-oo 



Gtj{t',t"\x) = Gc,{t' -t")+ / dtiGc.{t' ~h)U{t,ti\x)G1j{ti,t"\x). 



(2.36) 



Integrals over the entire Kadanoff-Baym-Keldysh contour survive only in the definition of Ga{t) 
in equation (|2.33p . which can instead be calculated directly on the real axis by picking t and t' on 
the upper branch of the contour in equation (|2.34p 



-t-oo 

G„(i'-t") = -^ J Img„(w)e-*"(*'-*") [/(w) - e(t' - i")] 



(2.37) 



Also, for further discussions we will need its Fourier transform 

4-00 



Ga(c^) 



J dte"^*Ga{t) = Rega(w) + »tanh^Img<^(w). 



(2.38) 



We are now ready to summarize the final formulas for the Green's function by evaluating all of 
the terms in equation (|2.3ip and then performing the integration over the parameter x. We begin 
with the formal matrix solution for Gij{t' ,t"\x) 



G^it',t"\x)^ di \Mtut2)-Ga.{h,t2)Uc{t,h\x)\\-' Go^itX) 



(2.39) 



Here Sc{t,t') is the Dirac delta function along the contour defined by J^dt'6c{t,t')f{t') — f{t). 
After substituting this result into equation (|2.27p and then integrating over the parameter x, we 
obtain for the greater Green's function the following results: 



(2.40) 



Gj{t)— — iwoexp< — i{Ef — fj,)t 



dx I dt' / di 



dlUt^t^jxl 

^LTo[t,t ) 

ft ax 



iwo exp < —i{Ef — ij)t — j dx Tr 



||I - GoUc(i|a:)|| Gq- 



dx 

— — iwfj exp {—i{Ef — fi)t + Indctc ||I — GoUc(t)||} 

= - zu.oe-''(^^^-^'* det, \\d,{h,h) - Go{ti,h)Uc{t, <2)|1 , 

where Tic denotes the trace of a continuous matrix operator expressed as a line integral over the 
contour and detc is the determinant of the continuous matrix operator defined on the contour. 
However, in our case, the time-dependent field Uc{t,t2) (and hence all nondiagonal components) 
are nonzero only on the upper real time branch of the contour for < ^2 < ^, which allows us to 
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reduce the problem to the determinant of a continuous matrix operator over the finite real time 
interval: 

G>{t) = -iwoe-'^^'f det[o^t] \\S{h - - Goih - t2)U\\ . (2.41) 

It is the restriction to this finite real-time interval that makes the numerical calculations much 
more efficient than other methods used previously. Note that all of the temperature dependence is 
in the Green's function Go, which is the time-ordered Green's function for the effective medium of 
the dynamical mean-field theory. 

In a similar fashion, we can derive the expression for the lesser Green's function {t = —t > 0): 

G<{t) = ™ie*(^^-^)*det[o,i] - t2) + Gi{ti - t2)U\\ . (2.42) 

Here, the Green's function Gi{ti — ^2) is defined in equation (j2.37p and is the time-ordered Green's 
function for the effective medium when there is an /-electron on the impurity site. 

When we perform numerical calculations to evaluate the continuous matrix determinants, we 
first discretize the time contour with a discretization step At, and then discretize the continuous 
matrix operator into a discrete matrix, which can be manipulated by standard computational 
libraries like LAPACK. The discretized representation for the Green's functions is then 

G>{t) = -^^«oe"'(^^"^^*det[o,t] ||% - W,Go{t, ~ t,)U\\ (2.43) 

and 

G<{t) = ™ie^(^/-^')* det[o,t] + W,Gi{t, ~ tj)U\\ , (2.44) 

where Wi are the quadrature weights for the discrete set of times {ti}; for a uniform grid Wi = 
W = At, we have to calculate the determinant of a Toeplitz matrix. This will allow us to use 
Szego's theorem and the Wiener-Hopf approach to find analytic expressions of the continuous 
matrix operator determinant which are accurate for long times. Note, that while it appears that 
these expressions are equally valid for all temperatures, the numerical solution, in the At — > 
limit, becomes more difficult at lower temperatures. 

The retarded Green's function is found from the greater and lesser Green functions via 

G}{t) = e{t){G>{t)-^G<{t)}. (2.45) 

Using the relation in equation (12. 8p . then yields 

G}{t) = -ze(t)e-'(^-f-^)* {wodet[o,t] ||I - Got/|| + (det[o,t] ||I + GiC/||)*} (2.46) 

for real times, and 

+ OC 

G}{co) ^-i j dte^('^+^-^/)* [wo det[o,t] ||I - GoC/|| + wi (det[o,t] ||I + GiL/H)*} (2.47) 


for real frequencies. 

3. The Wiener-Hopf sum equation approach and Szego's theorem 

We follow closely the development of the Wiener-Hopf sum equation approach in McCoy and 
Wu [ [20] , which is based on the original integral equation development by Wiener and Hopf [ [21] , 
as summarized in Krein [ I22j . The Wiener-Hopf integral equation approach has been used widely 
in the X-ray edge problem and in the Falicov-Kimball model [[S] primarily at T = 0; we will show 
that the Wiener-Hopf sum equation approach is very powerful for finite temperature calculations. 

Before we can apply this technique to the problem of calculating the /-electron spectral function, 
we must begin with some mathematical preliminaries (further details can be found in Ref. [[2U]). 
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Consider a sequence of numbers c„ which satisfy X]^-oo l*^"! ^ of + 1 hnear 

equations 

N 

m— 

for < n < N with c„ and j/„ known. We will eventually be taking the limit where ^ oo, so we 
will require J2'^=-oo I?/" I < too. Note that the matrix c„_m that appears in equation p.ip is a 
Toeplitz matrix, because it depends only on the difference n — m and not on n and m separately. 
Our job is to find the solution Xm ^ to these equations. We define 

xi""^ ^yn^O (3.2) 

for n < and n > N. In addition, we need to define two more sets of coefficients: 

„ - / J2Z=0 CN+n~mXm^ for n > 

" " ^ forn < 



and 



- / J2Z=0 C-n-mXm^ for 71 > 

for 71 < 



Next, we relax the condition on equation p.ip to apply it for all n and we allow the sum over tti 
to extend from — oo to oo (although all cases have just a finite number of terms because xiJ^^ is 
potentially nonzero for only N + 1 terms) . The result is 



+00 

E 

m— — 00 



W = y„ + e(„ > N)u„^N + e(n < 0)v-n- (3.3) 



We consider a variable ^ = exp[i6'] which lies on the unit circle, so that |^| = 1. We multiply both 
sides of equation p.3p by ^" and sum over n to find 



with 

+ 00 JV JV 

C(0= E C"^"' X(e)-5^a;Wr, nO-E2^"^"' (^.5) 

n— — 00 n— n=0 

00 00 

c^(o = E""^"' ^'^d no = E^"^"- 

n— 1 n—1 

It may appear that we have made the problem more complicated, because we have replaced a 
problem with one unknown X{£^) by a problem with three unknowns X(^), U{^) and V{^). But 
it turns out that this more complicated problem can actually be solved by carefully analyzing its 
analytic structure. 

The crucial step of Wiener and Hopf is to find a unique factorization of C(^) into two factors, 
one analytic inside the unit circle and continuous on the circle, and the other analytic outside 
the unit circle and continuous on the circle. Before we show how to do this, we need some more 
mathematical definitions. We call a function a + function if it can be expanded as a Laurent series 
for 1^1 — 1 (X)^o°'"'^") ^^"-^ coefficients satisfy J2'^=o ^ ^^^^ case, the function is 

analytic within the unit circle and continuous on the unit circle. Similarly, we call a function a — 
function if it can be expanded in a Laurent series for |^| = 1 (X]n=-oo the coefficients also 

satisfy X^ni-oo l*^"! This function is analytic outside the unit circle and continuous on the 

unit circle; it vanishes as |^| — > 00. Note that any function defined on the unit circle via a Laurent 
expansion /(f) = X]^-oo with J2'^=-oo I'*"! ^ °° "^^^ obviously be uniquely decomposed 
into the sum of a + and a — functions via /(C) = f+iO + f-iO- 
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The function C(^) is factorized in a so-called canonical factorization, where 

cio = p{0-'Q{r')-\ (3.6) 

and both P{£,) and Q{^) are nonzero for |^| < 1; the factorization is made unique by requiring 
(5(0) — 1. The functions P and Q can then be found in terms of + and — functions as 

P(^)=e^+(«) and Q(ri) = e«-(«\ (3.7) 

which automatically satisfies G_(oo) — [and Q{0) = 1] and makes both P(^) and Q{S,) nonvan- 
ishing inside the unit circle. An explicit calculation of G±(^) is nontrivial, but it is easy to write 
down a formal solution via 

G+(e)--[lnC(C)]+ and G_(e) = - [lnG(0]_ , (3.8) 

where the ± subscripts denote a restriction of the Laurent series for the logarithm to its nonnegative 
or negative power terms, respectively. Note that the condition for properly defining the + and — 
functions requires the logarithm to be single-valued after ^ winds around the unit circle, or, in 
other words, we require IndG(^) — 0, where 

IndG(0 = ^ [lnG(e2") - lnG(e*)] . (3.9) 

We will also need to consider situations where IndG(^) = —1. In this case, a new function C{^) = 
— ^G(^) will have zero index, and the canonical factorization can be applied to G(^) instead. 

We are now ready to solve equation (|3.4p by substituting in the factorization from equation p.6p . 
Next, multiply both sides of the equation by Q{£,~^) to get 

[pm-'x{o = Q{r')Y{o + Qir'wm" + Qic'wir'). (3.10) 

The term on the left hand side is a -I- function, because P{0) ^ 0, and [P(0]~^ can be expanded in a 
power series with just positive powers for |^| < 1 and X{£^) is a -I- function. Similarly, Qi£,~^)V{£_^^) 
is a — function. The other two terms are neither 4- nor — functions, but they can be uniquely 
decomposed into -I- and — function pieces. We do this, and move all -I- functions to the left hand 
side and all — function pieces to the right hand side. We are left with 

[Qir')Yio] _ + [QiOuioe] _ + Qic'wir')- (3.11) 

The left hand side of equation (|3.1ip is an analytic function for |^| < 1 and is continuous on the 
unit circle, the right hand side is an analytic function for |^| > 1 and is continuous on the unit 
circle, and both functions agree on the unit circle, due to the equality, so we can define a function 
that is analytic in the entire complex plane, and hence is an entire function. But this function 
vanishes as |^| — s- oo, and the only entire function that vanishes for large argument is the function 
that is identically equal to 0. Hence we learn that 

x{o = pio [Q{r')y{0] + + pio [Qir'wm''] + (3.12) 

and 

vic') - - [QiO] { [Q{r')Y{o] _ + [Qir'wio^]^} ■ (3.13) 

Using the fact that ^(C^^)C^ is a + function [recall X(^) is an order N polynomial in we 
can substitute ^ in equation (|3.1ip . multiply both sides of the equation by , and then 

separate into the + and — function pieces to create another vanishing entire function, and learn 
that both 

x(r^)^^ = Q(0 { [piOYir')^] + + [^^(r')^(c)c^] +} (3.14) 
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and 



(3.15) 



hold. These are all of the necessary relations we will need for determining the Toeplitz determinants. 

Now we must take another mathematical interlude to state Szego's theorem. This theorem 
determines how a Toeplitz determinant exponentially decays for large matrix size [j23j and even 
finds the constant prefactor for the exponential decay [ [211 US] . The proof can be carried out using 
the Wiener-Hopf sum equation approach [ [3D] , but it would take too much space for us to repeat 
it here, and it is rife with many technical mathematical manipulations that would take us away 
from the physical phenomena we wish to discuss here. So we will instead simply state the result, 
and then show how to evaluate asymptotic expressions for the /-electron Green's function in real 
time. 

If we consider the determinant 13 at of the N x N Toeplitz matrix defined by 



Co 


C-1 


C-2 


■ C-N+1 


Cl 


Co 


C-1 


■ C-N+2 


C2 


Cl 


Co 


■ C-N+3 


CN-2 


CN-3 


CN-4 ■ 


C-1 


CN-1 


CN-2 


CN-3 ■ 


Co 



then Szego's theorem says that 



lim D 



N- 



N 



exp 



Ngo 



^ ngn9-n 

n=l 



with 



dn 



1 

2^ 



^6*6-"^ lnC(e^''). 



(3.16) 



(3.17) 



The conditions that IndC(^) = 0, C(^) is continuous on the unit circle, and '}2n=-oo < oo are 
all required for the theorem to hold. 

Now, with all of the mathematical preliminaries complete, we are ready to apply these tech- 
niques to finding the Toeplitz determinants needed for the /-electron Green's function. We need 
to evaluate the two determinants in equations (|2.43p and (|2.44p . We will show explicitly how to 
calculate the determinant for the greater Green's function. Modifications for the lesser Green's 
function are straightforward. 

Case 1: no winding. We start by discretizing the time axis [0,t] with + 1 points given by 
tn = nAt, with At — t/N . Then we define coefficients c„ to satisfy 



c„ = 5no - AtUGoinAt). 



(3.18) 



Since Go{t) is bounded, and decays exponentially fast in time, it is easy to see that X]^J°=-oo l*^"! *^ 
oo. We will assume for the moment that Ind C(C) = (which turns out will be true at half filling 
when U < 0.866 on the hypercubic lattice, see below). Using the definition for c„, we find 



E 



c„e"^ = 1 - AtU J2 Go{nAt)t 



(3.19) 



If we let nAt — t' and ujAt — 9, we recognize that in the limit where At — s- 0, the sum in the right 
hand side of the above equation approaches the Fourier transform of Go{t'), so we write 



C{lu) = 1-U / dt'Go{t')e 



1-UGoito), 



(3.20) 
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where Go(w) is defined by equation (|2.38p . The coefficients 
written as 



in equation (|3.17p can then be 



git') = 



9n_ 
At 



1 

2^ 



T/At 



duie 



-inAtiU 



1 

2^ 



-TT/At 



where we again used t' = nAt. Now, defining Dpf to satisfy 

D{t) = D{NAt) ^Dn^ det[o,t] \\5,j - AtGo{U - tj)U\\ , 
we find from equation (|3.16p that 



lim D{t) — exp 

t—^oo 



^ J duj\nC{Lj) + j dt't'g{t')g{-t') 



(3.21) 



(3.22) 



(3.23) 



which is the asymptotically exact result in the limit of large t when the index of C(^) is equal 
to zero. Since the coefficient of t in the exponent can be a complex number, the determinant can 
oscillate while exponentially decaying at long times. This occurs when the system has undergone 
the Mott transition and the /-electron spectral function displays a gap at low frequency. 

But we can improve upon the result in equation (|3.23|) to provide finite-time corrections. We 
illustrate next how to do this using the Wiener-Hopf sum equation approach. We start with our 
original + 1 simultaneous equations we need to solve in equation (j3.ip , with the coefficients c„ 
defined in equation p.lSp . Next, we choose yn = ^„o- Cramers rule, applied to the first column of 
the c matrix immediately tells us that 



which further implies that 



Dn = 



n 

\M=N 



D 



N 



D 



N+l 



X Dm^oo- 



(3.24) 



(3.25) 



If we assume that x, 



(M) 




(1 — AISm) exp(— go) (which we will verify below), then we find 



D{t) = exp 



-At 6m 



M=N 





oo 






exp 


2^ / 






— OO 



dwlnC(cj)+ / dt't'g{t')g{-t') 



(3.26) 



This will then provide the next order of corrections to the asymptotic limit of the determinant. 

We need to find 2^0^^'' for large M. Our strategy is to set C/(^) — in equation (|3.13p to 
find T^(C~^)i then solve for C/(^~^) with equation p.lSp and the approximate T^(^), and finally 
substitute L/(^) into equation p.l2p to find X{S^). The term Xg*^'' then follows from taking the 
^ limit. Using the facts that Y{^) — 1 for our equation and (9(^~^) — 1 is a — function, we 
immediately find that equation (|3.13p can be solved by 



v{c') = -i+m-')Y' ■ 

Substituting into equation (j3.15p (after replacing ^ ^^^) then yields 



(3.27) 



(3.28) 



Now we must replace ^ — s- in the above equation. This requires care because we will convert 
the — function into a -I- function, but with no constant term in the power series expansion. So we 
write 



1 f-AT 



(3.29) 
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with the prime indicating there is no constant term in the + function. This result can now be 
substituted into equation (|3.12p to yield 



x{0 = PiO - PiO 



Qir'){Pim''c'' piO{Qir')} 'r^ 



(3.30) 



where we used the fact that [Q{£, ^)]+ = 1. We will need to evaluate this in the limit ^ ^ 0. 
Note that we have the following two identities for P and Q 



P{0 - exp 



and 



Q(r') = exp 



and P{lo) — exp 



and Q{—Lo) — exp 



dtg{t)e" 



-At 



dtg{t)e'' 



(3.31) 



(3.32) 



Using the definition for g{t) in equation (I3.2ip . one can immediately verify that P{uj)Q{—uj) = 
1/C{u!) in the limit as At 0, as it must. Using the first relation in equation (|3.3ip . immediately 
shows us that P{^ = 0) = exp[— go]- Evaluating the other term in equation p.30p requires more 
work. In order to create the relevant + functions, we must first take the functions of lu, Fourier 
transform them to functions of t and then reverse Fourier transform back, but include only the 
positive time contributions. For example, we can write 



pioiQioy'^ 



-N 



1 

2^ 



dcu 



At 



(3.33) 



The limit of At is to remind us that there should be no constant term in the expansion. Using 
this strategy, we can 
At(5jv) exp[— (7o] with 



this strategy, we can immediately write down the final answer for x"^^ . We obtain x'^^ = (1 — 



oo oo 



1 

4^ 



P{uj) Q(-w') 



(3.34) 



— C30 At — oo 

where t = NAt. Putting this all together finally yields the asymptotic expression for D{t) 



D{t) = exp 



OO OO 



^ j dwlnC(w)+ j dt't'g{t')g{-~t')^ j dt j dt' h{t + t')h' {-t - t') 

-oo t At 



with 



h{t) ^ — I duj 



P{^) 
-Qi-uj) 



and h'(t) = — duj 



P{u) ' 



, (3.35) 



(3.36) 



Case 2: winding of —1. Szego's theorem fails when the index of C(^) is nonzero, and InC(^) 
winds around the origin. This occurs on the hypercubic lattice at half filling when U > 0.866, see 
below. We can still derive asymptotic formulas for the Green's function by solving an auxiliary 
problem, which has the winding removed. We now show how this is done. 

In all numerical cases we have examined, the index of C(^) satisfies IndC(^) — —1 when U is 
large enough. Such a winding can be removed by considering C{^) = --^C(^) which has index zero 
(the minus sign is introduced for convenience, as will be clear below). The coefficients obviously 
satisfy 

c„ = -c„_i. (3.37) 
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We examine the auxiliary Toeplitz determinant 





Co 


C-1 


C-2 


■ C-N+1 




-C_l 


-C-2 


-C-3 


-C-N 




Cl 


Co 


C-1 


■ C-N+2 




-Co 


-C-1 


-C-2 


■ —C-N+1 


Dn = 


C2 


Cl 


Co 


■ C-N+3 




-Cl 


-Co 


-C-1 


■ —C-N+2 




CN-1 


CAr-2 


CN-3 ■ 


Co 




-CjV-2 


-CN-3 


-cn-a ■ 


-C-1 



and consider the + 1 simultaneous equations in equation 



with the barred variables 



N 



^-n—ni'^rn 



Vn 



(3.38) 



with yn = 6no- Then Cramers rule, evaluated for the N + 1st column, tells us 



N-ll 



(3.39) 



where the two factors of (—1)^ cancel. Szego's theorem can be applied to C(^), so we immediately 
learn that 



with 



lim Dn — exp 



1 

2^ 



Ngo 



OO 

n=l 



'lnC(e'''). 



(3.40) 



(3.41) 



We need to calculate a;^^ to find the determinant Dn- 

First note that equations (|3.12f[3T5|) hold for the barred functions when we have a winding 
around the origin. We will once again solve for t^(^~^) setting U{^) = 0, then substitute into 
equation p.l4p to find X{^~^)^^ . Taking the limit ^ ^ will then give us x^n\ Just like before, 
we find 

ViO = -l+[QiOr\ (3.42) 

which then gives 

xic')^'' = Qio [{Qioy' nr')^''] ^ ■ (3.43) 

We take the limit ^ — > by integrating the function on the right hand side over the unit circle and 
dividing by 2n. Using 6 — Atuj, then yields 



[-AT 



At 
2^ 



7r/At 



duje"^'P{-u})/Q{u}), 



(3.44) 



where we used t — NAt again. Note that the limits on the integration cannot be simply extended 
to ±oo as we did before. Since the function lnC(exp[i0]) winds once around the origin, it displays a 
discontinuity of — 2j7r to its imaginary part as 9 runs from — tt to tt. We remove this discontinuity by 
adding a linear function that is equal to at = — tt and is equal to 2j7r at 9 — t: (the minus sign in 
C is needed to move the branch cut of the logarithm from the positive real axis to the negative real 
axis) . Since this linear function cannot be extended to infinity, we need to work with a finite value 
of At in the calculations (for a further discussion, look at Fig. [2] and the corresponding discussion 
below) . We simply take At small enough, that we see the numerical results do not change, and the 
system has approached its At — s- limit. The final result for the determinant is then 



D{t) = exp 



t 

2^ 



dLj\nC{uj)+ / dt't'g{t')g{~t') 



At 
2^ 



Tr/At 



dw'e"^'*P{-u')IQ{uj'). (3.45) 



-ir/At 
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The limits on the first integral have been extended to ±00, because it turns out that the linear 
piece gives no contribution to the integral. The definitions of g{t), P{uj) and Q{uj) are obvious from 
the above discussion. 

It should be noted that, in contrast to the case without winding where the third term in the 
exponent of equation (|3.35p gives finite time corrections to the large t exponential asymptotics of 
equation p.23p . in the case with winding, equation (|3.45p has a time dependent prefactor with a 
nontrivial At and i — > 00 limit which can not be considered as a vanishing correction to the 
exponential decay of equation (|3.23p . 

4. Numerical results 

We begin our numerical discussion by describing the different behavior of C(^) when it has 
no winding or when it winds with an index of —1 around the origin. In Fig. [21 we show a plot of 
lnC{^) for two cases: (a) the case with no winding, and (b) the case with an index of —1. These 
two cases are distinguished by the sign of the real part of C(^) at 1^ = 0, where its imaginary part 
changes sign. On the hypercubic lattice at half filling we have ReC(O) > for U < 0.866 and there 
is no winding and ReC(O) < for [/ > 0.866 and now the index is equal to —1. Notice how the 
imaginary part has a steep change in its value near ^ = [panel (a)], which evolves into a discrete 
jump by 27ri when U > 0.866 [panel (b)], the jump at the origin is compensated by the linear 
shift which runs from the minimal to maximal values of the frequency used in the calculations 
(|a;| = IOtt here). While it is obvious that integration over to can easily be extended to u; = ±00 
when there is no winding, one needs to carefully include contributions present at the given value 
of At when there is winding, and one cannot extend the integration limits in this case. 



0.5 



3 

u 



-0.5 




I I I I I I I I I -JC I I I , I I 1 I I , L_ 

-20 -10 10 20 -30 -20 -10 10 20 

0) 0) 



Figure 2. (Color online) Plot of (a) InC(cj) for the case with no winding [U = 0.5, IndC(^) = 0] 
and (b) In C(lj) and InC'(cj) for the case with winding [U = 1.5, At = 0.1, IndC(0 = -1]. 

Next we compare our different asymptotic expansions to the results derived from a direct 
numerical calculation with the matrix formalism from equation (|2.43p when T = 0.1. This involves 
a straightforward approach where we use three different discretization sizes {At = 0.1, 0.0666, 
and 0.0333) which we extrapolate to At — > 0, and we compare those results to calculations with 
different approximations. In the case with no winding, we use both the asymptotic form for large t 
in equation p.23p and the more correct form for smaller t in equation (|3.35p . In Fig. [31 we show the 
results for the greater Green's function at half filling on the hypercubic lattice for U — 0.5 [panel 
(a)] and U = 0.7 [panel (b)]. Note how the exponential form, that comes from Szego's theorem, is 
quite accurate for large times, but becomes poor for small times, and how the corrections arising 
from the Wiener-Hopf approach produce remarkable agreement with the numerically exact results 
for essentially all t. The errors are the largest near t — 0, but even there, they lie below the few 
percent level. This shows that the analytic formulas are quite accurate for all t when there is no 
winding. 

Next we examine what happens in the case with winding. We examine two values of U in 
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Figure 3. (Color online) Real and imaginary parts of the greater Green's function for cases with 
no winding (a) U — 0.5 and (b) U — 0.7 (T — 0.1). The solid (black) lines represent exact 
results of the functional determinant in equation (|2.43|l . the dot-dashed (red) line represents 
the asymptotic result in equation (|3.23|l . and the dashed (blue) line represents the expression 
with finite-time corrections in equation (|3.35|l . Note that we cannot really discern any difference 
between the exact matrix results and the asymptotic Wiener-Hopf results all the way down to 
t = 0. 

Fig. m U = 1.5 which is near the critical interaction for the metal-insulator transition Uc = V^, 
and U ~ 2 which is a small gap insulator. We show the resuhs from the scaled matrix calculations 
with a solid line (black) and the asymptotic Wiener-Hopf results with a (blue) dashed line. In 
this case, the Green's function has behavior that does not approach the asymptotic exponential 
behavior for short times, and the analytic formula is less accurate for small times (with maximal 
error is on the order of 10%). 




Figure 4. (Color online) Real and imaginary parts of the greater Green's function for the case 
with winding [IndC(C) = -1] (a) U = 1.5 and (h) U = 2 {T = 0.1). The solid lines represent the 
extrapolated matrix results for the functional determinant in equation (|2.43|l and the dashed 
line represents the asymptotic expression with the finite-time corrections in equation (|3.45p . 



Now we can specify the differences between the cases without winding and with winding. 
For U < 0.866 there is no winding and our exact results of the functional determinants rapidly 
approach the asymptotic result (with finite-time corrections) in equation (|3.35p and we do not 
observe a crossing of the zero axis at any temperature and at T ^ the long-time behavior is 
replaced by a power law. On the other hand, for 0.866 < U < Uc we always observe a crossing of the 
axis at high temperatures, which originates from the time dependent prefactor in equation (|3.45p . 
With decreasing temperature, the crossing point shifts to larger times producing some kind of 
exponential decay for intermediate values of t and at T = we observe a crossover to power law 
behavior when the zero crossing is pushed out to infinity. For U > Uc the crossing of the axis is 
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observed at finite time values for any temperature and does not transform into power law at zero 
temperature. 

The Fourier transforms to the spectral formula produces results essentially equivalent to those 
shown already in Refs. [ UHl [H] 112] , so we do not repeat them here. 

5. Conclusions 

In this work, we have shown a new representation for the /-electron Green's function, which 
allows for efhcient numerical computation. By using the property that Uc{t,f) is nonzero only on 
the interval [0,t], we restrict the determinant of the continuous matrix operators to a finite time 
interval instead of over the three branches of the Kadanoff-Baym-Keldysh contour. This produces 
a huge savings in the computational effort, because the size of the matrices do not grow with 
temperature and they are less than one half the size of the matrices used in previous numerical 
calculations. As a result, we can examine much wider regions of parameter space. 

In addition, we used the Wiener-Hopf sum equation approach, along with Szego's theorem, 
to derive exact analytical expressions for the Toeplitz determinants required to find the Green's 
functions. While these expressions are asymptotically exact in the limit where t — > oo, we find they 
have errors less than 10% all the way down to t = 0. This then allows us to use the matrix results 
for small times, and then append them by the asymptotic expressions for large times in order to 
find accurate results for the Green's function at all times. These are then Fourier transformed to 
real frequencies to determine the /-electron spectral function. 

The /-electron spectral function displays the expected behavior as well. For U values smaller 
than the critical U for the Mott transition {U — \/2), the spectral function develops a sharp 
peak at low T, which ultimately diverges as an inverse power law in w as T — > [rigorously 
speaking, our asymptotic formulas are not valid at T = 0, because the function InC(^) becomes 
discontinuous due to the jump in the Fermi factor, but the power law behavior can be extracted 
via other techniques — our focus here was on the finite-T behavior]. When U is larger than the 
critical U, a gap forms in the /-electron spectral function (rigorously speaking, it is a pseudogap 
on the hypercubic lattice) but subgap states rapidly enter as a function of T for low T. 
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